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Abstract 

Anaerobic glycolysis in yeast perturbed by the reduction of xenobiotic ketones is 
studied numerically in two models which possess the same topology but different 
levels of complexity. By comparing both models' predictions for concentrations and 
fluxes as well as steady or oscillatory temporal behavior we answer the question 
what phenomena require what kind of minimum model abstraction. While mean 
concentrations and fluxes are predicted in agreement by both models we observe 
different domains of oscillatory behavior in parameter space. Generic properties of 
the glycolytic response to ketones are discussed. 



1 Introduction 



The reductionist attempt to identify the function of all genetic constituents of 
a living organism is one of the major goals of the post-genomic research [1,2]. 
In a living cell, however, single molecular constituents are linked together in 
a complex network including gene regulation, signaling and metabolic path- 
ways [3,4]. Whether the detailed knowledge of the individual properties of 
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metabolites taking part in enzymatic reactions may help to infer the large 
scale behavior of the living cell is another question challenging life-scientists, 
and its answer will necessarily imply a highly inter-disciplinary cooperative 
effort. 

A recent study of gene regulatory networks for early fruit fly development 
has suggested that the inherited robustness of the network's behavior may 
solely result from the network's topology [5]. This robustness may render the 
observed behavior insensitive to the quantitative details of the microscopic 
interactions. Hence modeling of such networks was successful even for rather 
crude approximations of the details. Other findings on artificial gene circuits 
showed a broad variety of different functions despite the same topology, sug- 
gesting that quantitative details become important for some topologies [6]. 
Moreover, the construction of a mathematical model and the corresponding 
abstraction of the complex intracellular network are always confronted with 
the problem of neglecting "unnecessary" details. The issue of the appropriate 
level of abstraction has been disputed particularly in the case of metabolic 
networks. Therein it has been suggested that modeling efforts risk to be no 
longer taken serious by experimentalists if modeling is continued by means of 
"too" strong abstractions or falsified toy models [7]. 

Here we compare predictions from two models of different levels of abstraction 
to characterize the consequences of neglecting details. For this fundamental 
study we choose the metabolic network of glycolysis (Fig. 1) in baker's yeast 
(Saccharomyces cerevisids) which is of particular interest to biotechnology [8,9] 
and models of different degree of complexity have been suggested and used 
[10,11,12,13,14,15,16]. We choose two models that translate the metabolic net- 
work into systems of ordinary differential equations (ODEs), a procedure that 
is well established [17,18,19]. Our subsequent analysis treats the ODEs fully 
nonlinear and focuses on the temporally asymptotic, self-organising behav- 
ior of the models. The distinction between steady and oscillatory states in 
the models appears to be crucial, since several physiological advantages of 
the latter have been proposed. The dependences of the concentrations and 
fluxes on external control parameters are numerically computed via a bifurca- 
tion analysis [20,21]. This approach enables detailed studies of large parameter 
spaces that would be more time consuming by simulation of the complete time 
courses including transient dynamics [22]. This issue will gain importance as 
the analysed theoretical models increase in size and complexity. 

In order to test the predictive power of both models, we perturb both mod- 
els by the same set of simple equations that account for the conversion of a 
xenobiotic (here we choose a ketone) by a single enzyme (here alcohol dehy- 
drogenase) into a chiral carbinol which consumes the cofactor nicotinamide 
adenine dinucleotide (NADH), see Fig. 1(a). This particular perturbation was 
chosen both for its pharmaceutical relevance and its broad impact on the back- 
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Fig. 1. Schematic view of a yeast cell in chemostat culture and models of different 
degree of abstraction, (a) Glucose influx at constant flow rate is parameterized 
by the concentration [Glc^Jo- The perturbation (dark shaded box) is controlled 
by the concentration [K x ] of ketone in the medium and described by the same 
set of equations with identical parameter values for two models. The full scale 
model (b) and core model (c) describe glycolysis (dotted box) at different degrees 
of complexity. Dots (lines) denote metabolites (reactions), "N" ("A") stands for 
NADH (ATP) and the light shaded box in (b) is considered in detail in Subsec. 3.2. 

bone of glycolysis via the NADH mediated interference with the redox balance 
of the cell. The pharmaceutical interest stems from the increasing demand 
for stereoisomerically pure pharmacologically active compounds (PACs). The 
stereo-geometry of PACs essentially determines the physiological effects of a 
chiral drug, as the thalidomide (contergan®) tragic showed in the 1950s [23]. 
Stereoisomerically pure carbinols are central precursors of modern, innovative 
pharmaceuticals and make for more than 50% of all essential substructures 
in pharmaceuticals [9]. Today stereo-selective carbinol synthesis by means of 
biotransformation, i.e. chemical reactions catalyzed by enzymes in vivo, is 
indispensable for the production of PACs [24,25,26]. Further progress in the 
understanding and rational manipulation of biotransformatory carbinol syn- 
theses strongly depends on mathematical models that account for the limiting 
reproduction of the cofactor (NADH in the considered case) and the underly- 
ing regulatory networks [27,28]. The models studied in this paper also consti- 
tute a first step in this direction. In the future we will extend the analysis to 
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other eucaryote's metabolic responses to a range of pharmacologically active 
compounds. 

We emphasize the exposure of systems in oscillatory states to xenobiotics as 
a sensitive experimental probe of the network's global response. The global 
response is experimentally accessible via NADH-ffuorescence, and oscillatory 
or steady behavior is easily distinguishable as individual cells synchronize their 
oscillations [11]. The temporal frequency of the oscillations can be measured 
with high accuracy. 

The paper is organized as follows. In Sec. 2 two models of glycolysis and 
the unique set of equations representing the perturbation are introduced and 
the employed numerical method of bifurcation analysis is outlined. Section 3 
presents our results which are discussed in Sec. 4. 



2 Models and Methods 

The glycolytic pathway of Saccharomyces cerevisice fermenting glucose to 
ethanol is one of the best studied metabolisms from both the experimental 
and modeling side [10,11,12,13,14,15,16]. We consider two different mathe- 
matical models in terms of ODEs that identically represent the topology of 
the glycolytic pathway under anaerobic conditions, see Fig. 1(b), (c). 

2. 1 Glycolysis 

Specifically we start from the full scale model of glycolysis devised by Hynne 
et al. [15] and compare its predictions with those of a core model with the 
same topology adapted from Wolf et al. [12]. Public domain "Silicon Cell" 
versions of the models can be found at http://www.jjj .bio.vu.nl, includ- 
ing a graphical overview and the mathematical expressions. The full scale 
model contains 22 variables for concentrations of involved metabolites and 
quantitatively accounts for most known details on enzyme regulation in order 
to precisely describe the supercritical onset of oscillations as observed experi- 
mentally [11,15]. The core model has been derived from the original 9 variable 
model [12] by addition of the variable Sf x (extracellular glucose) and the cor- 
responding flux balance of the medium = y yo ik([Glc x ]o — Sf x ) and a simple 
form of glucose transporter saturation j = k (Sf x — Si)/ (Sf° — S\ + K trans ) as 
well as lactonitrile formation jVcn = VvoikcTS(Sl x in the presence of cyanide. A 
storage flux jstore = ^storeSi^ is included in analogy with the full scale model 
but as a consequence of combining the three reactions downstream of glucose, 
istore is fed by intracellular glucose Si in the core model instead of glucose-6- 
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parameter 



value in present model value in model of Ref. [12] 



y,. n i ■ k ( min 1 ) 

yvoi v / 


59 • 0.048 


1.3 


fci (mM~ min ) 


6211.86 


100.0 


&2 (mM -1 min -1 ) 


38.3756 


6.0 


&3 (mM -1 min -1 ) 


126819 


16.0 


kr, (min -1 ) 


22.2857 


1.28 


k G (mM" 1 min -1 ) 


6.93908 


12.0 


k, (min -1 ) 


24.7 


13.0 


H 


4.0 


4.0 


Kj (mM) 


0.52 


0.52 


Jv total V imvi / 


0.98 


1.0 


A+n+ol (mM) 


3.6 


4.0 


fco (mM min -1 ) 


48.1 




-Ktrans (mM) 


0.002 




fcstoro (mM -1 min -1 ) 


17.1 




fccN (min -1 ) 


0.01477 




ADH total (mM) 


0.01 




fe^DH ( mM -2 min -l) 


80800 




feADH (min-l) 


1000 




k£ Dli (min -1 ) 


5400 





Table 1 

Comparison of parameters for our core model with the original model (Reference) 
by Wolf et al. [12]. The values of y vo i, k, n, N totsl , A total , ADH total , fe^ D 3 H are identical 
in both models and taken (except for ADH tota \ and fc^ D 3 H ) from Ref. [15]. For the 
full scale model, parameters are identical with Hynne et al. [15]. 



phosphate in the full scale model. The consumption of two subsequent entities 
ATP per glucose for storage would contribute with — 2j S torc = — 2/c Store S 1 74 3 
to the ATP balance. However, we find that this term overestimates the feed- 
back of a varying glucose concentration on the ATP balance and we always 
observed subcritical onsets of oscillations under this assumption. Indeed, not 
intracellular glucose but glucose-6-phosphate is the substrate and allosteric ac- 
tivator of the pathway to carbohydrate storage [29]. In the ATP balance, we 
therefore keep the ATP consumption for storage included in the term — k 5 A 3 
of unspecific ATP consumption in line with the original 9 variable model [12]. 
The remaining nomenclature and kinetic terms are identical with Ref. [12] but 
we have fitted a new set of parameters (see Table 1) to describe the same data 
as the full scale model. 

In both models, we have replaced the reduction of acetaldehyde (ACA) to 
ethanol (EtOH) by the same simple reaction scheme: 

k ADH 

ACA + NADH + ADH # 

k^ (1) 

£; ADH 

ADHaca 3 -> ADH + NAD+ + EtOH 



5 



which later enables the coupling with a second reaction that utilizes the same 
enzyme alcohol dehydrogenase (ADH; EC 1.1.1.1). We introduce the concen- 
tration [ADHaca] of the bound enzyme as an additional variable with the 
ODE 

d [ADHaca] _ + ( ) 

^ — y ADH,ACA y ADH,ACA ■ \ z ) 



Its gain (^adhaca) an d l° ss (^adhaca) fluxes are described by simple mass 
action kinetics in agreement with the Michaelis-Menten term used for the 
complete reaction in the original full scale model [15]. 

^adh,aca = kt DH ■ [ACA] • [NADH] ■ [ADH] (3) 
%dh,aca = A; 2 ADH • [ADHaca] + A; 3 ADH • [ADH AC a] (4) 

The unbound form of the enzyme is computed from the total amount ADH tota \ 
of the enzyme by 

[ADH] = ADH total - [ADHaca] (5) 

whereas the concentration of NAD + follows from 

[NAD+] = iV total - [NADH] - [ADH AC a] • (6) 



The ODEs for [NADH], [ACA] and [EtOH] have also been updated by the 
new expressions for the fluxes. We have fixed two of the four new parameters 
to ADH total = 10 fM and A; ADH = 1000 min and computed the remaining 
two self-consistently from the fluxes in the original full scale model [15], for all 
values see Table 1. This procedure preserved the full scale model's behavior 
as verified in Table 2 and by comparison of the bifurcation diagrams Fig. 2 (a) 
and Fig. 5 in Ref. [15]. Identical equations (2)-(6) and the same parameter 
values are used in both models. 

Although the core model possesses the same topology as the full scale model, it 
uses simple mass action kinetics for all lumped steps of the pathway [12]. One 
single regulatory interaction has been included, specifically the regulation of 
6-phosphofructokinase (PFK; EC 2.7.1.11) via ATP. This strong abstraction 
of the core model is sufficient to reproduce glycolytic oscillations [12]. 

We calculated a new set of parameters (see Table 1) for the core model such 
that the onset of oscillations at [Glc x ]o=18.5 mM is captured to the same ac- 
curacy (see Table 2) as in the full scale model, which had particularly been 
optimized to describe the onset of oscillations. To fit the core model, we in- 
serted the known values of fluxes and concentrations from the full scale model 
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Variable in core and full scale model Core model Full scale model Reference 



[Glc x ] 


[Glc x ] 


18.5 


18.501 


18.5 


gex 


[Glc x ] 


1.54982 


1.55409 


1.55307 


Si 


[Glc] 


0.560999 


0.573935 


0.573074 


S2 


[DHAP] 


2.08167 


2.97161 


2.95 


S3 


[BPG] 


0.000266672 


0.00026973 


0.00027 


St 


[ACA] 


1.48032 


1.49748 


1.48153 


cex 
D 4 


fACA 1 


1.28731 


1 30226 


1 28836 


N 2 


[NADH] 


0.330046 


0.329872 


0.33 


^3 


[ATP] 


2.08505 


2.09875 


2.1 


Jin 


jin 


48.0 


47.99 


48.0 


^storage 


^storage 


20.0 


19.90 


20.0 


jGlyc 


jGlyc 


4.77 


4.82 


4.8 


jEtOH 


jEtOH 


46.46 


46.54 


46.4 


joutACA 


joutACA 


4.77 


4.82 


4.8 



Table 2 

Comparison of metabolite concentrations (top, all values are in mM and accurate 
to six significant digits) and fluxes (bottom, rounded and in mM min -1 ) for both 
models at the onset of oscillations. Reference values are taken from a comprehensive 
modeling effort that incorporated available experimental data (Table 6 in [15]) and 
only a subset of the full scale model's variables are shown. All others match the 
Reference equally well. 

(Table 6 in Ref. [15], here included as Reference in Table 2) into the expres- 
sions of the core model and calculated those parameters kj which are factors 
of mass action terms. This strategy was partly used in Ref. [15]. Parameters 
ki of flux terms that involve a yet undetermined saturation constant were 
expressed as functions (with inserted values of fluxes and concentrations at 
onset) of the respective saturation constant. Then a steady state at the onset 
of oscillations (Hopf bifurcation) was determined for an initial choice of the 
unknown parameters and numerically continued along curves in the space of 
unknown parameters (see Methods below). Note, this strategy is very efficient 
since each run explores two unknown parameters simultaneously which are 
linked by the constraint of reproducing a Hopf bifurcation. We monitored the 
values of the resulting fluxes and fixed that combination of remaining param- 
eters with closest correspondence of the computed fluxes to the data from the 
full scale model (Tables 7,8 in Ref. [15]). The computed concentrations at the 
onset of oscillations in the core model were also in good agreement with those 
of the full scale model (see Table 2). The temporal period of the oscillations 
near onset was close to 0.65 min in the full scale model and 0.25 min in the 
core model. 

Fig. 2 shows that the new set of parameters also yields a supercritical onset of 
the oscillations in the core model in agreement with the full scale model and 
experiments [11,15]. This set of parameters is therefore well suited for further 
studies of synchronisation phenomena in a multicellular context [13]. Note, 
that both models also share a second oscillatory domain at lower glucose con- 
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[GlcJ in mM [GlcJ in mM 

Fig. 2. Bifurcation diagrams of (a) the full scale model and (b) the core model of 
glycolysis showing steady states (thin curves) as well as minimum and maximum 
(thick) of NADH oscillations for varied inflow glucose concentration [Glc^Jo- Solid 
(dashed) curves represent stable (unstable) states. 

centration. This has also been observed in the original full scale model (Fig. 8 
in Ref. [15]) but so far not in experiments. This agreement of both models' 
predictions may stem from their identical topology but still is remarkable since 
it covers a parameter range outside of that where the models had been fitted. 

Altogether, both the full scale model and the core model show very similar 
behavior and create the impression that quantitative details are largely irrel- 
evant when the model is analysed in the same context that had already been 
used to determine the remaining unknown parameters. Hence "predictions" 
from a simplified model are reliable for conditions similar to those considered 
in the model's construction, e.g. both models predict a second oscillatory do- 
main for low glucose supply. However, we want to test the reliability of the 
core model under conditions that may be very different from those included 
during parameter optimization. Such reliability is essential for the rational 
improvement of biotechnological processes beyond the conditions applied so 
far. 

2.2 Perturbation of the redox-balance 

Both the full scale model and the core model have been augmented with 
the same set of equations that account for the biotransformation of a xeno- 
biotic ketone. In the considered case the ketone ethyl acetoacetate is re- 
duced by alcohol dehydrogenase (ADH; EC 1.1.1.1) to enantiomerically pure 
ethyl L-3-hydroxybutyrate (carbinol). Although this reaction has been stud- 
ied extensively, little is known about how the intracellular enzymes inter- 
act [30,31,32,33,34,35]. Currently, two enzymes, D-directing /3-ketoacyl reduc- 
tase (KAR; EC 1.1.1.100) and L-directing ADH are known to be involved in 
the asymmetric reduction of the substrate [28,36]. Yet, KAR activity is negli- 
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gible with a large excess of sugar present. Hence, the overall stereo-selectivity 
of the bioconversion largely depends on cell physiology [28]. 

In analogy to the reduction of acetaldehyde (ACA) in reaction (1), we now 
consider the reduction of ketone (K) to carbinol (C) by the same simple reac- 
tion scheme: 



K + NADH + ADH 



J, ADH 
ft 6 



i,ADH 

"<4 



(7) 



ADH K ADH + NAD+ + C 



with rate constants k^ Dn , k^ Dn and /c^ DH . If not stated differently then the 
new parameters are chosen equal to the corresponding values of the acetalde- 
hyde reduction, i.e. A^ DH = dh ^adh = ^.adh and ^adh = ^adh The 

corresponding ODE with gain (f adh k) an d l° ss ( v adh k) h uxes reads 



d[ADH K ] _ + _ 

^ — ^AD^K — ^ADH.K \°) 

« W = kt DU ■ [K] • [NADH] • [ADH] (9) 
«W = ^5 ADH " [ADH K ] + A; 6 ADH • [ADH K ] . (10) 

The additional variable [ADH K ] also alters the constraints and yields via the 
consumption of [NADH] a strong feedback on the upstream steps of glycolysis. 



[ADH] = ADH total - [ADHaca] - [ADH K ] (11) 
[NAD+] = iV total - [NADH] - [ADH AC a] - [ADH K ] (12) 

The membrane transporters for ketone Vk = k ■ ([K x ] — [K]) and carbinol 
were chosen linear in order not to restrict the accessible range of intracellular 
perturbations. The permeability for ketone k = 24.7 min" 1 was set equal to 
the value for acetaldehyde [15]. 

The above perturbation introduces a competition between acetaldehyde and 
ketone for unbound enzyme ADH and the cofactor NADH. In Sect. 3, the 
behavior of both models is tested for consequences of the perturbation and 
predictions on oscillatory or steady behavior are compared. 
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2.3 Methods 



The investigated metabolic network may be considered as a nonlinear dynam- 
ical system that evolves its own state, i. e. the set of values for the metabolites' 
concentrations, in a self-organizing manner. Various computational methods 
have been developed and proven to be essential for the analysis of complex 
cellular networks [18,19], e.g. bifurcation analysis [20,21] and direct numerical 
simulation [22]. 

Since we focus on steady and oscillatory states it is most convenient to combine 
the computation of the system's state along "branches" in order to efficiently 
scan the parameter space. This strategy is called "continuation" and starts 
from a single known solution and proceeds in subsequent steps to compute 
similar solutions for neighboring values of the parameters, e.g. [Glcjo (mixed 
flow glucose concentration at inlet), [K x ] (extracellular ketone concentration). 
We used the software package AUTO for such tasks [37] . Therein the solutions 
are represented as the root of an extended system of algebraic equations de- 
rived by discretizing the ODEs. Any new solution (the root) is then computed 
iteratively by the Newton algorithm. 

The software also monitors the linear stability of the solutions and thereby 
detects critical values of the parameters where the solutions change drastically, 
so called "bifurcations". The onset of oscillations, a Hopf bifurcation, is of 
special interest here hence the employed strategy amounts to a bifurcation 
analysis. Bifurcations may themselves be continued in two control parameters 
simultaneously which was used to fit the parameter set of the core model to 
flux and concentration data of the full scale model at the Hopf bifurcation. We 
have supplemented the AUTO package by subroutines for oscillatory states 
that monitor and export the time averages of fluxes and the relative phase 
shifts between periodic time courses of concentrations. 



3 Results 

3. 1 Comparison of model predictions 

The fluxes within the metabolic network of glycolysis are redirected when 
the extracellular ketone concentration is increased. Both models show almost 
identical behavior, see Fig. 3. The effluxes (storage, ethanol, glycerol, acetalde- 
hyde) are represented by the intervals between the curves above and below the 
corresponding area and all effluxes add up to the total influx in units of C3. The 
observed increase in the efflux of acetaldehyde corresponds to the additional 
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Fig. 3. Flux diagrams (log-linear scales) at fixed glucose concentration 
[Glca;]o=30mM for (a) the full scale model and (b) the core model. All fluxes have 
been measured in equivalents of C3 entities and have been normalized by the in- 
flux (96.2 mMcj min" 1 ). Individual pathways (storage, ethanol production, glycerol 
production, acetaldehyde (ACA) efflux as annotated) consume a portion of the in- 
flux given by the distance between the curves above and below the respective area. 
(c),(d) The biotransformation flux for both models. 



consumption of cofactor NADH via ketone reduction. The latter is fostered by 
increasing ketone concentration, and the stoichiometric constraint on NADH 
recycling (light shaded area equals sum of dark shaded areas) demands an 
equal flux of NADH regeneration. Glycerol secretion is discriminated by the 
lack of NADH, however, in the full scale model a small flux to glycerol remains 
even for large ketone concentration whereas this flux is strongly suppressed in 
the core model. We observed the same behavior independently of the value of 
[Glcjo- Fluxes have been averaged over time for the oscillatory states. Note, 
the time courses of substrates and cofactors in oscillatory states are generally 
out of phase which influences the average fluxes, see Subsec. 3.2. 

The dependence of the oscillations on the strength of the ketone reduction 
is shown in Fig. 4. Both models behave differently with sustained oscillations 
predicted by the full scale model over the whole range of perturbations whereas 
the core model does favor steady behavior. However, the temporally averaged 
concentration of NADH shows the same depression upon ketone increase, re- 
flecting the additional NADH consumption. 

We repeated the analysis for several values of the carbinol release rate /c^ DH 
from the ADH complex, see Eq. (7). Different release rates are expected for 
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Fig. 4. Bifurcation diagrams for (a),(c) the full scale model and (b),(d) the core 
model at fixed [Glc x ]o=30mM. Oscillations of the NADH concentration (gray area in 
(a),(b)) occur between a maximum and a minimum value (solid curves) that decrease 
for stronger perturbation by ketone. (c),(d) Temporal period of the oscillations for 
k£ DR = 5400 min 1 (solid curve) that decreases (frequency increases) with [K x ] for 
the full scale model in (c). Additionally, data for fcj^ DH = 4800 min" 1 (dashed) and 
fcg DH = 3800 min" 1 (dotted) indicate that periods may also remain unaffected or 
increase with [KJ depending on the carbinol release rate. 



structurally dissimilar ketones due to varying interactions between the enzyme 
and structures besides the common keto-group. The corresponding variability 
in the affinity k^ DH and in A;^ DH are not considered as such changes can be 
compensated by rescaling of [K x ] and do not yield new results. 



Fig. 4 (c) shows how the oscillation period is affected by ketone for three 
choices of the release rate /c^ DH . The oscillations speed up (period decreases) 
for /c^ DH > k^ u = 4800 min -1 including the case of equal release rates 
fcg- DH = kg Dn and the oscillations slow down (period increases) if the release 
rate is chosen slower than k^ n . The critical rate /c^? H arises from the compe- 
tition of additional NAD + recycling flux due to additional substrate (ketone) 
versus diminished NAD + recycling via ACA reduction when the enzyme is 
trapped longer in the ADH K complex and therefore less active. The former 
effect accelerates the oscillations whereas the latter slows them down. Sus- 
tained oscillations have been observed over long time scales [11] which enable 
accurate measurements of the oscillation period. Such measurements can be 
used to infer rate constants of different ketones relative to the rates of ACA 
reduction and to distinguish different ketones. 
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Fig. 5. Phase diagram of stable stationary (white) and oscillatory (gray) states in 
the [Glc x ]o - [K x ] parameter space (linear-logarithmic scales) for (a) the full scale 
model and (b) the core model. 



Fig. 4 (d) confirms that the core model does not possess the dynamical fea- 
tures discussed above. In order to characterize the conflicting predictions of 
both models, we have computed the parameter values of the onset of oscilla- 
tions in a two-parameter plane, see Fig. 5. We observe a qualitative difference 
in terms of the relations between the three Hopf bifurcations that are present 
in the unperturbed system. The full scale model possesses extended oscilla- 
tory domains (gray shaded). When [K x ] is increased, two Hopf bifurcations 
at large [Glc x ] are connected and the enclosed steady behavior vanishes. The 
core model favors stationary states: two Hopf bifurcations at small [Glc x ] are 
connected and the enclosed oscillations vanish when [K x ] is increased. Hence, 
the different levels of abstraction mainly determine the dynamical aspects of 
the models' behavior. On the other hand the averaged or stationary quantities 
are much less affected by the considered abstraction. 

We further examined the case of low membrane permeability for ketone or 
equivalently low external concentration of ketone (to maintain viability) which 
both yield small intracellular concentrations of ketone. At a fixed ketone con- 
centration ([K x ] < 1 mM), we varied [Glc x ]o and monitored the fluxes. A 
maximum of carbinol production is found at an intermediate value of [Glc x ]o 
which depends on the fixed value [K x ] (not shown). Both models predict this 
maximum which can be interpreted as follows. The lack of NADH recycling is 
limiting the carbinol flux for small [Glc x ] whereas at large [Glc x ] almost all 
ADH binds to acetaldehyde which is then much more abundant than the com- 
peting ketone. Hence at intermediate value an optimum carbinol flux develops 
where both limiting influences are weak. In biotechnological applications it is 
often desired to adjust the glucose supply accordingly. 
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3.2 Phase shift in oscillatory states 



Finally, we have closer analysed possible differences between temporally av- 
eraged oscillatory fluxes and the coexisting steady state fluxes, both in the 
full scale model. Note, the steady states are unstable when oscillations occur 
supercritically but unstable states can nevertheless be computed by the con- 
tinuation algorithm as this method does not rely on temporal integration. The 
unstable steady state represents a solution that equally respects stoichiomet- 
ric and regulatory constraints as the stable oscillatory state does and a minor 
change in parameters can render the steady state stable. Hence the coexisting 
unstable steady state reproduces the behavior of an alternative operational 
mode. 

As an example, we compare the fluxes towards glycerol production for a range 
of values of the carbinol release rate /q^ DH from the ADH complex, see the 
remarks above and Eq. (7). Fig. 6 (a) shows a small portion of the reaction 
scheme of glycolysis where dihydroxyacetone phosphate (DHAP) and NADH 
are reaction partners. The concentrations of fructose 1,6-bisphosphate (FBP), 
glyceraldehyde 3-phosphate (GAP) and DHAP perform synchronized oscil- 
lations whereas the cofactors NADH and NAD + always possess anti-cyclic 
oscillations. Therefore, this branch point of glycolysis (Fig. 6 (a) and shaded 
area in Fig. 1 (b)) may act as a switch if the oscillation of either NADH or 
NAD + can be synchronized with the substrate pool. 

The computations are performed in the full scale model and all parameters 
of the glycerol path are kept fixed. Fig. 6 (b) shows the averaged oscillatory 
glycerol flux that is increased by more than 50% at large rate /c^ DH with 
respect to the steady state at the same rate. The time courses of substrate 
DHAP and cofactor NADH are depicted in Fig. 6 (c,d) for two different values 
of /c^ DH . The temporally averaged concentrations (c) t are almost identical 
with the steady state value at the same /c^ DH . We define a phase 0(c) = 
27r(t(c max ) + t{c m in)) /2T to characterize a time course with period T and 
time of maximum (t(c max )) and minimum concentration (t(c m j n )). The large 
increase in averaged oscillatory flux versus the steady state value stems from 
the oscillations of DHAP and NADH that are almost in phase at large /c^ DH 
(Fig. 6 (d,e)). 

Slower release of NAD + (smaller kg Dn ) increases the phase shift (A0 = 
0([DHAP]) - 0([NADH])) to the trailing DHAP oscillation (Fig. 6 (c)) and 
decreases the averaged oscillatory flux relative to the steady state flux (Fig. 6 
(e)). In addition to the phase shift A0, the averaged oscillatory flux is also 
affected by changes (due to varied /q^ DH ) in amplitude and functional form of 
the oscillations. Moreover, the absolute values of both steady and averaged os- 
cillatory flux vary strongly (Fig. 6 (b)) as the temporally averaged metabolite 
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Fig. 6. (a) Small portion of the reaction scheme of glycolysis corresponding to the 
shaded box in Fig. 1(b). The substrate DHAP is reduced by NADH to finally yield 
glycerol, (b) Temporally averaged glycerol flux v osc in oscillatory state (thick curve 
between two Hopf bifurcations at k$ DU ps 2250 min -1 and 4750 min -1 ) is found to 
be larger or smaller, respectively, than the corresponding steady state flux Steady 
(thin curve), (c) Oscillatory time courses of [NADH] (thin) and [DHAP] (thick 
curve) are out of phase by Acp = 1.2 at /c^ DH =2300 min -1 . Average concentrations 
(c)t have been subtracted and [NADH] has been scaled by a factor of 10 in (b,c). (d) 
Both concentrations oscillate almost in phase (A(f> = 0.2) at larger fcj^ DH , here 4600 
min -1 . (e) The relative difference of fluxes (v osc — v s teady)/^steady (thick) is partly 
due to varying phase shift (thin curve) between the time courses of substrate 
DHAP and cofactor NADH. Parameters are [Glc x ] = 15mM and [K x ]=10mM. 

concentrations change with /c^ DH . The biotransformation flux towards carbinol 
is affected by changes of similar absolute size which, however, amount to slight 
relative changes only. 



In the left half of the oscillatory interval, the averaged oscillatory flux to 
glycerol is even smaller than in steady state which is only possible for a large 
phase shift, as observed in Fig. 6 (c). The effect of the phase shift becomes 
dominant for large oscillation amplitudes because it clearly matters, if the 
reaction partners frequently collide when they appear at the same time or if 
they only rarely react when one of them is always too late. This situation is 
analogous to a problem in electrical engineering where the electrical power at 
a device in a sinusoidal alternating current (AC) depends on the amplitudes 
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of voltage and current and on the phase shift between their maxima. 



4 Conclusions and perspective 

We have compared predictions of two theoretical models that have previously 
been used to describe glycolytic oscillations in yeast. Both models possess the 
same topology but incorporate different levels of kinetic details. Irrespective of 
these differences, both models predicted identical dependencies of temporally 
averaged fluxes and concentrations on the strength of the perturbation, i. e. the 
extracellular concentration of a xenobiotic ketone. This observation provides 
confidence in theoretical modeling given the topology of the network is known 
and completely incorporated. Note, the unknown parameters of the models 
were fitted to data of the unperturbed system hence the models predicted 
the response to the perturbation. The chosen perturbation affected the entire 
network via multiple feedback loops of cofactors. We suggest the exposure to 
xenobiotics as a feasible experimental probe of the network's global response. 
In the future we will extend this analysis to metabolic responses in other 
eucaryotes and to a range of pharmacologically active compounds. 

The kinetic details of both models affected the temporal organisation of the 
responses. Oscillations ceased to exist in the abstract core model while the 
full scale model supported oscillations for a large range of perturbations. Os- 
cillatory behavior has recently been suggested by Hauser et al. to protect 
enzymes and possibly cells against toxic substances [38,39]. Oscillations also 
play an important role for calcium signal transduction [40] where they ap- 
pear in bursts [41,42]. Bursting, the temporally repeated sequence of quies- 
cent and oscillating episodes, has recently been found in a theoretical study 
of stochastically driven glycolysis [43] . Oscillations have also been suggested 
to improve the thermodynamic efficiency of glycolysis with respect to steady 
operation [44]. Furthermore, kinetic details are important for (non periodic) 
short term responses to glucose pulses [10]. 

A closer inspection of the flux distribution in oscillatory states revealed that 
temporally averaged fluxes may differ from the values in the coexisting (un- 
stable) steady state, albeit for each metabolite the (unstable) steady state 
concentration equals the temporally averaged value in the oscillatory state. 
This effect results from the relative phase shift between time courses of sub- 
strates and cofactors particularly at branch points of the network. The time 
average of any oscillatory flux is generally restricted to a convex combination of 
elementary flux modes [45]. However, allosterically regulated enzymes impose 
additional restrictions on the fluxes at given metabolite concentrations which 
can result in suboptimal convex combinations particularly in steady states. 
Oscillatory states possess additional degrees of freedom, i. e. phase shifts, that 
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are adjusted by many even remote reactions. Evolution of metabolic networks 
may therefore have acted faster on those networks capable of oscillatory be- 
havior enabling them to explore more possibilities and to become superior. 
This conjecture is in line with the observed oscillatory behavior in one of the 
most important and inherited metabolic pathways, i.e. glycolysis, and it sup- 
ports the view that oscillatory behavior may be more abundant than presently 
recognised [38,39]. 



A corresponding strategy for metabolic engineering may target the phase shifts 
in oscillatory states as has been suggested by Ross and coworkers [46]. They 
forced simple reaction schemes with oscillating substrate concentrations to 
increase desired fluxes. Here we found that the same effect is also important for 
large metabolic networks where phase shifts are adjusted in a self-organizing 
fashion. The artificial addition of appropriate buffers or bypasses (as in the 
present example) in order to control phase shifts are promising candidates for 
metabolic engineering. 



In terms of modeling strategies, one should treat those approaches with caution 
that focus on steady states, e.g. canonical modeling [47], since the computed 
flux distribution in an unstable steady state does not necessarily approxi- 
mate the actual oscillatory behavior around that steady state. Continuation 
algorithms as applied in the present paper appear to be particularly suited 
to address the (dis) advantages of oscillatory versus (hypothetical) stationary 
behavior. 



To summarize, we have shown that static features of the metabolic network of 
glycolysis are robust to crude abstractions of models but more care has to be 
taken when dynamic features are considered. Many open questions remain and 
call for further research work. E.g. the possibility to conduct ketone-to-carbinol 
biotransformations under aerobic conditions could increase the NADH recy- 
cling and hence the carbinol production which, per se, is an important biotech- 
nological challenge. Further modeling will help to solve the related problems 
by simulating the effect of suppressing or enhancing parts of the complex 
metabolic network. 



We are indebted to Markus Bar, Sune Dan0, Jean-Marie Frangois, Ursula 
Kummer, Karl-Heinz van Pee, Preben Graae S0rensen and Jana Wolf for fruit- 
ful discussions. MB wishes to thank Gerhard Jorg for practical assistance and 
the Graduate Research Center " Heterocycles" as well as the German National 
Merit Foundation for financial support. 
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